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Abstract 

In this paper, we are interested in the numerical approximation of the classical time-dependent drift-diffusion system 
near quasi-neutrality. We consider a fully implicit in time and finite volume in space scheme, where the convection- 
diffusion fluxes are approximated by Scharfetter-Gummel fluxes. We establish that all the a priori estimates needed 
to prove the convergence of the scheme does not depend on the Debye length A. This proves that the scheme is 
asymptotic preserving in the quasi-neutral limit A — * 0. 



1. Introduction 

1.1. Aim of the paper 

In the modeling of plasmas or semiconductor devices, there is a hierarchy of different models: kinetic models and 
quasi hydrodynamic models, ranging from Euler-Poisson system to drift-diffusion systems (see 1 38, 39[ 32]). In each 



of these models scaled parameters are involved, like the mass of electrons, the relaxation time or the rescaled Debye 



length. There is a wide literature on the theoretical validation of the hierarchy of models (see |4),|33j,ll2|] and references 
therein). Moreover, an active and recent field of research consists in designing numerical schemes for these physical 
models which are valid for all range of scaled parameters, and especially when these parameters may tend to 0. These 
schemes are said to be asymptotic preserving. These methods have proved their efficiency in many situations, for 
instance: in fluid limits for the Vlasov equation, quasi-neutral limits for the drift-diffusion, Euler or Vlasov equations 
coupled to the Poisson equation, in diffusive limit for radiative transfer (see HEISTS Hill [Tot] among a long 
list of articles that could be mentioned here) 

In this paper, we consider the numerical approximation of the linear drift-diffusion system. It is a coupled system 
of parabolic and elliptic equations involving only one dimensionless parameter: A, the rescaled Debye length. This 
parameter A is given by the ratio of the Debye length to the size of the domain; it measures the typical scale of 
electric interactions in the semiconductor. Many different numerical methods have been already developed for the 
approximation of the drift-diffusion system; see for instance the mixed exponential fitting schemes proposed in 10] 
and extended in B3lLl3411 to the case of nonlinear diffusion. The convergence of some finite volume schemes has been 
proved by C. Chainais-Hillairet, J.-G. Liu and Y.-J. Peng in JHIa]. But, up to our knowledge, all the schemes are 
studied in the case A — 1 and the behavior when A tends to has not yet been studied. 

In this paper, we are interested in designing and studying a scheme for the drift-diffusion system applicable for any 
value of A. This scheme must converge for any value of A > and must remain stable at the quasi-neutral limit A — * 0. 
We consider an implicit in time and finite volume in space scheme with a Scharfetter-Gummel approximation of the 
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convection-diffusion fluxes. As it is classical in the finite volume framework (see [18]), the proof of convergence of 
the scheme is based on some a priori estimates which yield the compactness of the sequence of approximate solutions. 
In the case of the drift-diffusion system, the a priori estimates needed for the proof of convergence are L°° estimates on 
N and P, discrete L 2 (0, T, //'(Q))-estimates on N, P and 4* in the non-degenerate case |Ht], with additional weak-BV 
estimates on N and P in the degenerate case J9fl. However, the crucial point in our work is to establish that all the a 
priori estimates do not depend of A > and therefore the strategy used in [8, 9] to get them does not directly apply. In 
order to get estimates which are independent of A, we adapt to the discrete level the entropy method proposed by A. 



Jiingel and Y.-J. Peng in [33] and by I. Gasser et al in [24, 25]. The choice of the Scharfetter-Gummel fluxes for the 



discretization of the convection-diffusion fluxes is essential at this step. 
1.2. The drift-diffusion system 

Let Q be an open bounded subset of R d (d > 1) describing the geometry of a semiconductor device and T > 0. The 
unknowns of the linear drift-diffusion system are the density of electrons and holes, N and P, and the electrostatic 
potential W. In this paper, we assume that the doping profile vanishes in the semiconductor device. The system writes 
for all f) e Qx [0,T]: 

d.N + div(- VN + 2VW) = 0, (la) 
d,P + divi-VP - PVY) = 0, (lb) 

- Pay = p - n, (lc) 



where A > is the rescaled Debye length. The system is supplemented with mixed boundary conditions (see [38]): 
Dirichlet boundary conditions on the ohmic contacts and homogeneous boundary conditions on the insulated boundary 
segments. It means that the boundary dQ. is split into dQ = T D UF N with r^nr*' = and that the boundary conditions 
write: 

N(y, t) = N D {y), P(y, t) = P D (y), ¥(y, t) = W D (y), (y, t) e T D x [0, T], (2a) 
(VN ■ v) (y, t) = (VP ■ y) (y, t) = (VP • v) (y, = 0, (y, t) e T N x [0, T] , (2b) 

where v is the unit normal to dQ outward to £1 

The system (Q]) is also supplemented with initial conditions A^o, Po- 

N(x, 0) = A^o(jc), P(x, 0) = P (x), x e O. (3) 

In the sequel, we denote by (Pj) the drift-diffusion system (|TJ— 0J- We need the following assumptions: 

Hypotheses 1.1. The domain Q is an open bounded subset ofW 1 (d > I) and dQ. — T D U Y N with Y D n T N — and 
m(r D ) > 0. The boundary conditions N D , P D and ^V are the traces of some functions defined on the whole domain 
Q, still denoted by N D , P D and Y . Furthermore, we assume that 

M),P(>eL°°(Q), (4a) 

N D , P D e L°° n H\Q), e H l (Q), (4b) 
3m > 0, M > such that m < N , P , N D , P D < M a.e. on O. (4c) 

The weak solution of CP A ) is defined by: N, P e L°°(Q x (0, T)), N - N D , P - P^, T - e L°°(0, T\V), with 
V = {v e H l (D.) ; v = almost everywhere on T D } and, foraU test functions <p e C~(Qx[0, T)) and 77 e C t °°(Qx(0, T)) 
such that <p(y, t) = i](y, t) = for all (y, f) e Y D x [0, T) and: 

f f (Nd t <p-VN-V<p + NV¥-V<p)dxdt+ f N (x) <p(x, 0) dx = 0, (5a) 

(Pd t (p- T VP-V(p-PTP¥-V<p)djcdt+ P (x)<p(x,Q)dx = 0, (5b) 

Jo Jn ^ 2 f |vF.Vj7djcA= f f (P-N)r]dxdt. (5c) 
Jo Jn Jo Jn 

The existence of a weak solution to the drift-diffusion system (P^) has been proved in [ 4(J ^]] under hypotheses more 



restrictive than Hypotheses 11.1 I since they consider more regular boundary conditions. In 02311 . the authors prove this 
existence results under Hypotheses 1 1.11 and assuming that V(logA r o - *Po), V(logPfl + *Po) are in L ro (Q). 
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1.3. The quasi-neutral limit of the drift-diffusion system 

The quasi-neutral limit plays an important role in many physical situations like sheath problems 1I20I1 . plasma diode 
modeling [44], semiconductors fl45ll ... Then, it has been studied for different models: see II 1 31. 14311 for the Euler-Poisson 



model, [6, 28] for the Vlasov-Poisson model and [33, 24, 25] for the drift diffusion-Poisson model 



In these models, the quasi-neutral limit consist in letting tend to the scaled Debye length A. In a physic point of view, 
this means that only the large scale structures with respect to the Debye length are then taken into account. Formally, 
this quasi-neutral limit is obtained by setting A = in the model, here (Pf). Then, the Poisson equation (ITcb on *F 
reduces to the algebraic relation P - N — (which is the quasi-neutrality relation). But adding and subtracting dTab 
and fib) , we get new equations on N and T*. The quasi-neutral system (Po) rewrites finally for all (x, t) e Q x [0, T]: 

d,N - AN = 0, (6a) 
div(jVVT0 = 0, (6b) 
P = N. (6c) 

In i33ll . A. Jiingel and Y.-J. Peng performed rigorously the quasi-neutral limit for the drift-diffusion system with a zero 
doping profile and mixed Dirichlet and homogeneous Neumann boundary conditions. Indeed, under Hvpotheses ll.il 
and under quasi-neutrality assumptions on the initial and boundary conditions (No - Pq = and N D - P D = 0), they 
prove that a weak solution to (Pf), denoted by (N A , P A , v P /l ), converges, when A -> 0, to (N°, P°, T* ) solution to (Po) 
in the following sense: 

N A -> N°, P A -> P° in Z/(Q x (0, T)) strongly, for all p s [1, +00), 

N A N°, P A -> P°, *P A -> T* in L 2 (0, T, H X (CIJ) weakly. 

The same kind of result is established for the drift-diffusion system with homogeneous Neumann boundary conditions 
by I. Gasser in ll24ll for a zero doping profile and by I. Gasser, C.D. Levermore, P. Markowich, C. Schmeiser in 1E5I1 
for a regular doping profile. In all these papers, the rigorous proof of the quasi-neutral limit is based on an entropy 
method. 

The entropy method, described for instance in the review paper (11], has been developed in the last twenty years. It 
is firstly devoted to the study of the long time behavior of some partial differential equations or systems of partial 
differential equations and to the study of their equilibrium state. It consists in looking for a nonnegative Lyapunov 
functional, called entropy, and its nonnegative dissipation, connected within an entropy-entropy production estimate. 
Generally, it provides the convergence in relative entropy of the evolutive solution towards an equilibrium state. This 
method has been widely applied to many different systems: see [1] and the references therein, but also [[35, 22, HI 27, 



26] 



However, the entropy method also permits to get new a priori estimates on systems of partial differential equations 



via a bound on the entropy production, see 0331 l24l 12511 for instance. In the case of Dirichlet-Neumann boundary 



conditions, the entropy functional, which has the physical meaning of a free energy, is defined (see 03311 ) by 
E(f) = j (h(N) - H(N D ) - \og(N D )(N - N D ) 

A 2 \ 
+H(P) - H(P D ) - \og(P D )(P - P D ) + — I VP - VT' D | 2 j(/x, 

with H(x) = j. log(f) dt — x log x - x + 1 , and the entropy production functional is defined by 

1(0 = J (^V|V(log^V- 1 I')| 2 + P|V(logP + >I')| 2 ) dxdt. 

The entropy-entropy production inequality writes: 

dE 1 

-r(t)+-Kt)<K D Vr>0, (7) 
dt 2 

where Kq is a constant depending only on data. This inequality is crucial in order to perform rigorously the quasi- 
neutral limit. Indeed, if E(0) is uniformly bounded in A, (0 provides a uniform bound on J I(s)ds. It implies a priori 
uniform bounds on (A^ 1 , P A , v F i ) solution to (pf) and therefore compactness of a sequence of solutions. 
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1.4. Presentation of the numerical method 

In order to introduce the numerical scheme for the drift-diffusion system (Pj), fi rst , we define the mesh of the domain 
Q.. Here, we consider the two-dimensional case but generalization to higher dimensions is straightforward. The mesh 
M = (T, £, P) is given by T, a family of open polygonal control volumes, £, a family of edges and P = {xk)k£T 
a family of points. As it is classical in the finite volume discretization of elliptic or parabolic equations with a two- 
points flux approximations, we assume that the mesh is admissible in the sense of 11 80 (Definition 9.1). It implies that 
the straight line between two neighboring centers of cell (xk, x£) is orthogonal to the edge cr = K\L (and therefore 
collinear to vk,o-, the unit normal to cr outward to K). 

We distinguish in £ the interior edges, cr = K\L, from the exterior edges, cr c <9Q. Therefore £ is split into £ = 
U &exi- Within the exterior edges, we distinguish the edges included in Y D from the edges included in T N : 
£>ext = &ext u &exf For a gi ven control volume K e T, we define &k the set of its edges, which is also split into 
£>k = £>K,im U £% ext U £>x exr For all edges cr e £, we define do- = d(x K ,x L ) if cr = K\L e £■„, and d^ = &(x K ,cr) if 
cr e & ext with cr s & K . Then, the transmissibility coefficient is defined by Tq- = m(cr)/do-, for all cr s £. We assume 
that the mesh satisfies the following regularity constraint: 

3£ > such that d(x K , cr) > £ diam(tf), VK eT,Vcre8 K . (8) 

Let p > be such that card(£/j-) < (3. Let At > be the time step. We set N T = E(T I At) and f = nAt for all 
< n < Nt- The size of the mesh is defined by size (T) = max^ £ 7- diam (K) with diam(/T) = sup A . yeK \x - y\, for all 
K e 7~. We denote by 6 = max(Af, size (T)) the size of the space-time discretization. Per se, a finite volume scheme 
for a conservation law with unknown u provides a vector uj- = (uk)k£T e (with 6 = Card(T)) of approximate 
values and the associate piecewise constant function, still denoted u-j-: 

where \k denotes the characteristic function of the cell K. However, since there are Dirichlet boundary conditions 
on a part of the boundary, we need to define approximate values for u at the corresponding boundary edges: u^d = 
(Uo-)o-eS d 6 (with 9 D = Card(£f r ,)). Therefore, the vector containing the approximate values in the control 
volumes and the approximate values at the boundary edges is denoted by um = 0*r, wgo). 
For any vector um = (ur, «£ X we define, for all K e 7", for all cr e Sk, 

if cr = K\L € G KJnt , 

ifo-6 6* (10a) 




and Da-u - YDukA ■ (10b) 
We also define the discrete H l - semi-norm | • on the set of approximations by 



O-Gfi 



As we deal in this paper with a space-time system of equations (Pa), we define at each time step, < n < Nt, the 
approximate solution u n T = {u" k )ket f° r u = N, P, "P and the approximate values at the boundary = (mJ-)^^ 
(which in fact does not depend on n since the boundary data do not depend on time). Now, let us present the scheme 
that will be studied in the sequel. First, we discretize the initial and the boundary conditions. We set 

(N°,F%)=-^f (N Oc),Po(d)dx, VK € T, (11) 

(jv* f* £(^(y), P°(r), ^ D (r)Vr, Vcr e £f„. 

fJ = P2, n=*£. Vcre£f AY ,Vn>0. (12) 



and we define 
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This means that N n & , = Ng D for all n > 0. 

We consider a Euler implicit in time and finite volume in space discretization. The scheme writes: 

m(K) K — K + J] T~Zt = °> V R e T ' Vn ^ °> ( 13a ) 

creE K 

pn+\ _ pn 

m(K) K — K + Y &"k)t = 0. v ^ e 7". v " ^ °' ( 13b ) 

-^Yj T °- DV K,<r = ™(K)(P" K - N n K ), SK e T, V« > 0. (13c) 

o-e£ s 

It remains to define the numerical fluxes T'^ and Q"^ which can be seen respectively as numerical approximations 

of I (- VN + JVW) • Vk,(t and J (-VP-PVT*) • v K ,<r on the interval [f", f" +I ). We choose to discretize simultaneously 

the diffusive part and the convective part of the fluxes, by using the Scharfetter-Gummel fluxes. For all K e T, for all 
cr e & K , we set: 

= T - (B(-£>^)^ +I - B(D^)^) , (14a) 

where B is the Bernoulli function defined by 

5(0) = 1 and B{x) = — - Vx + 0. (15) 

exp(jc) - 1 



These fluxes have been introduced by A. M. II' in in [29] and D. L. Scharfetter and H. K. Gummel in 14211 for the 
numerical approximation of convection-diffusion terms with linear diffusion. It has been established by R. Lazarov, 
I. Mishev and R Vassilevsky in 13611 that they are second-order accurate in space. Moreover, they preserve steady- 
states . With this choice of numerical fluxes, M. Chatard in 111 111 proved a discrete entropy estimate, with control of 
the entropy production, which yields the long-time behavior of the Scharfetter-Gummel scheme for the drift-diffusion 



system. The generalization of these fluxes to nonlinear diffusion has been studied by A. Jiingel and R Pietra in [34], 
R. Eymard, J. Fuhrmann and K. Gartner in lfl7ll and M. Bessemoulin-Chatard in fl. 

Remark 1.2. Let us note that the definition d 1 0b ensures that D^"'^ = and also that T'^ — Q' 1 ^ — 0, for all 
o~ G £>K ext - These relations are consistent with the Neumann boundary conditions d2bl >. 

In the sequel, we denote by (S£) the scheme (TTTt — dT3Tl . It is a fully implicit in time scheme: the numerical solution 
(N'x +l , P"% 1 , x ¥"k 1 )k£T at each time step is defined as a solution of the nonlinear system of equations <TT3T> — (fl~4T >. When 
choosing instead of D'V"^ in the definition of the fluxes (TBI) , we would get a decoupled scheme whose solution 

is obtained by solving successively three linear systems of equations for N, P and *P. However, this other choice of 
time discretization used in Qlgt] induces a stability condition of the form At < CA 2 (see for instance 0[). Therefore, 
it cannot be used in practice for small values of A and it does not preserve the quasi-neutral limit. 
Setting A = in the scheme (Sa) leads to the scheme (So). The scheme for the Poisson equation ( 1 1 3ct becomes 
P" K - N" K = for all K e T, n e N. In order to avoid any incompatibility condition at n — (which would correspond 
to an initial layer), we assume that the initial conditions No and Po satisfy the quasi-neutrality assumption: 

Po-M)=0. (16) 

Adding and subtracting (II 3ab and (I13bl >. we get 

K A( K + ^ Z + &&) = °' V * e T > V " * °> 

<teS k 

and J] - ffg) = 0, VK e T, Vn > 0. 

o-e£ a - 
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But, using P" K = N" K for all K e T and neN, and the following property of the Bernoulli function 

B(x) - B(-x) = -x VxeR, (17) 

we have, € T, Vcr € & KJ „, U & N Kext : 

and + g£ = -tv (B(DV n £) + B(-Z)«P#)) D1V£ . 

Let us note that these equalities still hold for each Dirichlet boundary edge cr e £^ e „ if N® = P®. In the sequel, when 
studying the scheme at the quasi-neutral limit (So), we assume the quasi-neutrality of the initial conditions ( TToT l and 
of the boundary conditions: 

P D -N D = 0. (18) 
Finally, the scheme (So) can be rewritten: VK €.T,Vn> 0, 

N" K +l -N n K ^ B(DV£i)+B(-Dy£i) 
m(K) K - K - £ ^ —DN n £ = 0, (19a) 

- J] TrBVgW? 1 + AO = 0, (19b) 

cteSk 

P n K -N n K = Q, (19c) 
with the initial conditions (TTTt and the boundary conditions (fl2ll , 

7.5. Mam results and outline of the paper 

The scheme (Sa) is implicit in time. Then we begin by proving that the nonlinear system of equations cTT~3T > admits a 
solution at each time step. The proof of this result is based on the application of Brouwer's fixed point theorem. The 
existence result is given in Theorem ll.3l and is proved in Section|2] 

Theorem 1.3 (Existence of a solution to the numerical scheme). 

We assume Hypotheses 17.11 let T be an admissible mesh ofQ satisfying © and At > 0. If A = 0, we further assume 
the quasi-neutrality of the initial and boundary conditions (1161 > and (118l >. Then, for all A > 0, there exists a solution 
to the scheme (Sa)-' (N'^, P^^'^ket e (R 9 ) 3 for all n > 0. Moreover, the approximate densities satisfy the following 
L°° estimate: 

V7C e T, Vn > 0, m < N n K , P\ < M. (20) 

Then, in Section [3] we prove the discrete counterpart of the entropy-dissipation inequality (|7). As the functions N D , 
P D , X V D are given on the whole domain, we can set: 

«, p° K , ?»°)= ^ x^ D(x) ' pD(x) ' ^^y*' s,k e r - 

For all n e N, the discrete entropy functional is defined by: 
E" = J] m(K) (H(N" K ) - H(N§) - log(A^) (N n K - N%)) 

KeT 

+ £ MK) {H(P\) - H(P D K ) - \og(P D K )(P\ - P D K )) + t. \¥ n M - V D M \[ M , 

KeT 

and the discrete entropy production is defined by 

«e* + min [P" K , P n Ka ) [Do- (log 7>" + y S n ) ) 

The discrete counterpart of (|7) is given in Theorem 1 1.4 1 



Theorem 1.4 (Discrete entropy-dissipation inequality). 

We assume Hypotheses 17.71 let T be an admissible mesh of £2 satisfying ([8} and At > 0. Then, there exists K E , 
depending only on Q., T, m, M, N D , P D , ^i 10 , f3 and ^ such that, for all A > 0, a solution to the scheme (S,i), 
(N!j-, P^-, ( P^-)o<«<A' 7 -, satisfies the following inequality: 

- F" 1 

— + - < K E , Vn > 0. (21a) 

At 2 

Furthermore, ifN° and P° satisfy the quasi-neutrality assumption (116t . we have 

N T -1 

^ Afl" +1 < K E (l +A 2 ). (21b) 

»=o 



Let us note that the last inequality ( 121bt . which ensures the control of the discrete entropy production, depends on 
A. However, as we are interested in the quasi-neutral limit A — > 0, we can assume that A stays in a bounded interval 
[0, A max ] and then get a uniform bound in A. 

In Sectional we show how to obtain, from the discrete entropy-dissipation inequality, all the a priori estimates needed 
for the convergence of the scheme. These estimates are given in the following Theorem 1 1.5 1 There are weak-BV 
inequality (12251 and L 2 (0, T, T/^-estimates tG^ on TV and P and L 2 (0, T, //^-estimates (|22cT > on 



Theorem 1.5 (A priori estimates satisfied by the approximate solution). We assume Hypotheses \l.l\ let T be an 
admissible mesh o/Q satisfying (O and Af > 0. We also assume that the initial and boundary conditions satisfy the 
quasi-neutrality relations d!61 > and (118b . Then, there exists a constant K E depending only on Q, T, m, M, N D , P D , ^f , 
P and % such that, for all A> 0, a solution to the scheme (S^), {N 1 ^-, P" T , v F^-)o<«<A' r , satisfies the following inequalities: 



ii=0 cre£ n=0 o-efi 

N T -\ 



J] AtY J T<TD W n+1 kD ( rP n+l f + (ZVV" +1 ) 2 )< K F (1 + A 2 ), (22a) 

ii=() creE 

-1 N T -1 

Yj A( Tj t -t( d ctN" +1 ) 2 + At Z T ^ P " +l f ^ *H1 + (22b) 

n=0 o-efi 

Vt--1 

J At J] T a {D^" +l ) 2 <K F (l+ A 2 ). (22c) 

n=0 o-efi 

Estimates (122bl i and i22c\ yield the compactness of a sequence of approximate solutions, as shown for instance in fl, 
applying some arguments developed in II 180 . To prove the convergence of the numerical method, it remains to pass to 
the limit in the scheme and by this way prove that the limit of the sequence of approximate solutions is a weak solution 
to {Sa)- It can still be done as in [8], but dealing with the Scharfetter-Gummel fluxes as in (5J. The convergence proof 
is not detailed in this paper. Let us just note that the convergence proof holds for all A > 0. 

Finally, in Section [5] we present some numerical experiments. They illustrate the stability of the scheme when A 
varies and goes to 0. They show that the proposed scheme is an asymptotic -preserving scheme in the quasi-neutral 
limit since the scheme order in space and time is preserved uniformly in the limit. 



2. Existence of a solution to the numerical scheme 

In this Section, we prove Theorem 11.31 (existence of a solution to the numerical scheme (Sf) for all A > 0). As 
the boundary conditions are explicitly defined by (fT2t . it consists in proving at each time step the existence of 
(7V^-, P n T , *P^-) solution to the nonlinear system of equations ( fT3l when A > or (T% when A — 0. We distinguish 
the two cases in the proof. 
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2.1. Study of the case A > 

We consider here A > 0. The proof of Theorem 1 1.3 l is done by induction on n. The vectors N^- and i*. are defined by 
(fTTT l and by (1 1 3cl >. Furthermore, the hypothesis on the initial data (|4cT > ensures that 



m < iV° , P° K < M \/K eT. 



We suppose that, for some n > 0, (A^, P n T , "Pi.) is known and satisfies the U° estimate ( 1201 ). We want to establish the 
existence of (A 7 ^ 1 , P^ 1 , x ¥p' ) solution to the nonlinear system of equations (Q~3), also satisfying d2"0l . Therefore, we 
follow some ideas developed by A. Prohl and M. Schmuck in 114 111 and used by C. Bataillon et al in ||2[]. 
Let n > 0, we introduce an application T£ : I 9 x I 9 -> M 9 x K e , such that T^(N T ,P T ) = (N T ,P T ), based on a 
linearization of the scheme ( TT3l and defined in two steps. 

• First, we define ¥7- € R e as the solution to the following linear system: 

-A 2 Yj TaDV^ = m(K)(P K - N K ), VK e T, (23a) 
with ¥o- = V<x e 6f«. (23b) 

• Then, we construct (Ay, P7-) as the solution to the following linear scheme: 

+ J] T<r (B (-D'JV) Nk - B (DV^) N Ko .) = 0, Vtf e 7", (24a) 

+ J] r (T (fi (OT^) Pj, - B /v) = 0, V^e T, (24b) 

with N<r = N° and P a = P D a Vcr e fif„. (24c) 

The existence and uniqueness of *iy solution to the linear system (1231 are obvious. The second step d24"l > also leads 
to two decoupled linear systems which can be written under a matricial form: Ay = S'L and ApP-j- = §" p . The 
matrix for instance is the sparse matrix defined by 

(A*r)*r = ^ (l + + Yj T <rB(-DV K/r ) VKeT, 

(Aivk.L = -TV B (OTir i0 -) VL e T such that <x = K\L e £,■„,. 

We verify that A^ has positive diagonal terms, nonpositive offdiagonal terms and is strictly diagonally dominant with 
respect to its columns. It implies that A# is a M-matrix: it is invertible and its inverse has only nonnegative coefficients. 
The same result holds for Ap. Thus, we obtain that the scheme d24"l > admits a unique solution (Ay, Pt) £ R e x R e , so 
that the application T£ is well defined and is a continuous application. 

Now, in order to apply Brouwer's fixed point theorem, we want to prove that !Tjj preserves the set 

C m , M = {(AV, P r ) eWxW; m< N K , P K < M, VK e 7"} . (25) 
The right hand side of the linear system (124a) is defined by 

<re£&.„ 
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If N T > 0, we have S" N >0 and, as A^ is a M-matrix, we get N T > 0. Similarly, if P r > 0, we obtain that P T > 0. 
In order to prove that Nk < M for all K e T, we introduce M7- the constant vector of R e with unique value M and we 
compute An(Nt - M7-). Using the property (1171 1. we get that for all K e T, 

1 — \ m(X) „ m u 

(A w (AT r - Mr)), = -^-Wl - M ) + ^ § (A ^ ~ M) 

Since fi is a nonnegative function and N D satisfies d4cb . we have, for all <x e £^ , 

B (Z)^V) 2\£ - B (-ZW^) M - B (DW^) (N° - M) - M < -DV^ M 

Then, using the induction assumption N" K < M for all K e T, it yields 

(A*(/? r - M r )), < ™^ ^ft-M)-M^r ff OTV, 

and using ( 1231 ). we get that for all A" G 7~ 

/ ~ x m(/Q / u \ m(/Q 

{A N (N T - Mr)), < - MJ (TV* - M) + M -^(fi - M )- (26) 

We can prove exactly in the same way that, for all K e T, 

(A N (N T - m T )) K > (£ - m) (JVjc - m) + m ^t^k - m). (27) 

Since // > is an arbitrary constant, we can choose it such that mAt < M At < fi. Then, if (Nt,Pt) e C,„,m, 
inequalities d2oT ) and d27l > imply that 

A N (N T - M r ) < and AjvCNr - n*r) > 0. 

As A w is an M-matrix, we conclude that m < N K < M for all K eT. The proof that m < P K < M for all 7v e 7~ is 
similar and we have {Nt,Pt) e Cm,M- 

Finally, r" is a continuous application which stabilizes the set C m ,M- Then, by the Brouwer's fixed-point theorem, T" 
has a fixed point in C m> M which is denoted by (A^- +1 , P^ 1 ) and satisfies the L°° estimate (l20b . The corresponding 4V 
defined by g3) is denoted by v P^ fl and C^ 1 ,^ 1 ,^ 1 ) is a solution to the scheme O- This shows Theorem [T31 
when /I > 0. 



- M Yj TvDl!^ 



2.2. Study of the case A — 

Now, we prove Theorem 11.31 when A = 0. In this case, thanks to the quasi-neutrality assumptions, we have shown 
that the scheme {So) rewrites as the nonlinear system of equations (TT9V Indeed, it is sufficient to study the system 
(fT9ab-(fT9bT>. whose unknowns are (A^ +1 , »F^ +1 ). 

The proof is done by induction as in the case A > 0. Let us first note that Nfj- satisfy the L°° estimate ( l20b . Then, 
we assume that, for n > 0, N^- is known and also satisfies (l20b . As in the case A > 0, we introduce an application 
T" : (R* ) e -> M fl such that T n {N T ) = N T , based on a linearization of (I19ab - (l 19bb and defined in two steps. 

• First, we define TV e K e as the solution to the linear system: 

- J] t^^ANk + Nk,o-) = 0, VA" e T, (28a) 

cr£E, K 

with »Fv = 1 ¥° V<x e «5^„. (28b) 
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• Then, we define Nr s K 6 as the solution to the linear system: 

m (K) - „ vn BiD^Ko-) + Bi-D^Ka) - 

-^■W* -N" K ) - J] t„-± K ' a ' 2 ^DNk^ = 0,VKe T, (29a) 

with A^ = A£ Vcr e & D ext . (29b) 

First, let us prove that the application T" is well defined. If Nk > for all K e T, the matrix of the linear system (128} 
is a positive symmetric definite matrix (it can be proved for instance by multiplying (12 8 at by and summing over 
K e T). Therefore, *¥j- is uniquely defined. 

The linear system d29l can be written under the matricial form A^AV = S" N where the matrix A# is defined by: 
(A^^^P + i J] T(r (B(D^ K ^) + B(-D^ K , a )) VK e T, 

o-e£s-\£^ 

(Ajv)^ = -y (BtOT^) + B(-D^)) VLef such that cr = K\L e £,„,. 
and the right hand side S^, is defined by: 

The matrix A/y is an M-matrix because it has positive diagonal terms, nonpositive off diagonal terms and it is strictly 
diagonally dominant with respect to its columns. Therefore the linear system d29l has a unique solution AV> so that 
the application T" is well defined. It is also continuous. 
Now, let us prove that T" preserves the set 

'KmM = {n t e K e ; m < N K < M, VK e 7~} . 
Therefore, we compute An(N-j — M7-). We obtain 

(A N (N r - M r )) K = " AO 

+ \ J] r (r (B(OT^ (r ) + B(-OT^ (r ))(<-AO. 



Thanks to the induction hypothesis and (Hell , we deduce that An(N<-/ — M7-) < 0. Similarly, we prove that An(N^, — 
ni7-) > 0. This implies AV £ 'Km.M- We conclude the proof of Theorem 1 1.3 1 in the case A — by applying the 
Brouwer's fixed point theorem as in the case A > 0. 



3. Discrete entropy-dissipation inequality 



In this Section, we prove Theorem 1 1.41 Therefore, we adapt the proof done by M. Chatard in [11] for the study of the 
long-time behavior of the scheme (in this case, the entropy functional is defined relatively to the thermal equilibrium). 
Since H is a convex function, we have E" > and E" +1 - E" < T\ + T2 + T3, with 



T\ = J]m(^)(log(A^r I )-log(^))(A r r 1 -^)^ 

KeT 

T 2 = J]m(K)(log(l^ 1 )-log(p^)(P^ 1 -Pl), 



KeT 



T, - — |vt/«+1 _ U)D I 2 _ ^_ Imn _ mD I 2 
1 3 ~ 2 I M m u,m 2 ' M M " 
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Multiplying the scheme on N ( 1 1 3ab by At(log\N'^ 1 ) - log [N^jj, summing over K e T and following a discrete 
integration by parts (using (fT2li). we rewrite T\ : 

r, = -At J] 2 r£! (log (ivf 1 ) - log (<)) 

= A* J] ((D log N" +1 )k,o- - (D log JV D w) . (30) 

cteE 
<teE k 

Starting from the scheme on P d 1 3ab and following the same kind of computations, we also rewrite T?. 

T 2 = At J] &" K ^((DlogP" + ') K ^-(DlogP D ) K ^). (31) 

o-eE 
(teE k 

Now, in order to estimate T3, we subtract two consecutive time steps of the scheme on *P ( 1 1 3cl >. It yields: 
-A 2 J] T„B¥g + A 2 Y, ^ = ™(K)((P'k ' " n) - ' - 

Thanks to the schemes on N d 1 3ab and P d!3bK it rewrites 

A 2 J] r,(D^ - DWlJ -A 2 Y, t,A.D9^ - DS*,) = A* £ i&& ~ 

cteEk cteEk o~eEk 

Multiplying this equality by ^t"^ 1 - X P^, summing over K eT, integrating by parts and using the boundary conditions, 
we obtain: 



A 2 J] TjpVg - D^ D Ka f -A 2 J] rADr' Ka - DV D K<T )(DV" K ^ - D^ D K(T ) = 



a-eE 
cteEk 



cteE 
a-eEa 



a-eE 
oeE k 



But, for all K e T and all <x e Sk, we have 



and therefore for all A > 



o-eE 
o-efi. 



From (00), (ED and (03, we get 



A/ < J] !F£! (D(log JV - Vfg - DQogN - 

creE L 
o-eE k 

+ ffg (D(log P + - DQog P + »P)y 



(32) 



But, thanks to inequalities ( IA.3al ) and dA.3bb . we have 



z 

creE 
o-eSk 
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Now, using ( IA.4al ). ( IA.4bb and Young's inequality, we get 



2 



2 max(A^ +1 ,A^) 2 



mm(Nl +l ,N'^)\D tT (logN - + 



mm(N" K +1 ,N n K ^) 



T -|D 4r (logJV-'P) i 



max^ 1 ,/^) 2 



+ min(P^\ Fjg) \DA\og P + ¥)" 
Finally, thanks to the L°°-estimates (l20b in Theorem ll.3l we obtain 



min(/» +1 ,J»£) 



HT> - E" _ 1 , , D D , 2 , D D{2 

? 9m I™ * ^M\i,m 2m 1 M Ml1 



which rewrites 



Af "2 

E n+1 _ E « , 



2m 
M 2 



Ar 



(33) 



But, thanks to hypothesis ( Bbb . the functions log(W D ) - *P D and log(f D ) + m D belong to //'(^)- Therefore, using 
Lemma 9.4 in 01811 . we have 

|log(< ) - «P^£ „ < 7C || log(JV D ) - V D f Hl 



(O) 

and |log(P° ) + < 7C || log(P D ) + >F% 



(£21 



with 7C depending on /? and £ (defined in ([8])). It yields ( 121 at 
Summing ( 121 at over n€|0,...Wr-l) yields: 



JVr-l 



J] Afl" +1 < + ^ Afl" +1 < + E°. 



(34) 



17 = 



n=0 



It remains now to bound E°. As the function H satisfies the following inequality: 

1 {y-xf 



\fx, y > H (y) - H(x) - log x(y - x) < 



we get, using d4ct . 



min(x, y) 2 

(M - m) 2 



X MK) (H(N°) - H(N%) - \og(N°) (< - N°)) < m(O) 



2m 



and the same inequality for P. Then, multiplying the scheme ( 1 1 3cb at « = by T*^. - *P2 and summing over K e T, 
we get 

^ 2 J] T lT DV K JDy° KiT - D^J = ^ m(*D(< - iV°)CP° - T°) = 0, 



cre£ 



ATer 



if the initial conditions satisfy the quasi-neutrality assumption ( TToT l. Then, using a(a — b) > (a — b) 2 /2 - b 2 /2 for 
a,b eR and once more Lemma 9.4 in II 1811 . we obtain 



,1 



A 



Y^m'^mUm < yl^l < Y^ll^llfliffi,. 
with 7C depending on /3 and £ (defined in (|8]l). 

Finally, we obtain E° < K®(1 + A 2 ), with depending on Q., m, M, T , f3 and Inserting this result in (l34l . we 
deduce the discrete control of the entropy production (12 lbt with an adaptation of the constant 

It concludes the proof of Theorem 1 1.41 Let us note that the hypothesis on the vanishing doping profile is not directly 
necessary to follow the computations in this proof. However, we need it in order to ensure the lower and the upper 
bounds on the discrete densities, with their strict positivity. 
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4. A priori estimates on the scheme 



This Section is devoted to the proof of Theorem ll.5l This proof is splitted into three steps: first, we establish the weak- 
BV inequality on N and P (122al >; then, we deduce the L 2 (0, T, //^-estimate on N and P (I22bb : finally, we conclude 
with the L 2 (0, T, ^-estimate on ¥ d22cT >. 



4.1. Weak BV-inequality on N and P 

First, let us first prove the inequality (122at of Theorem 1 1.51 Therefore, we denote by Tbv the left-hand-side of ( 122al >. 
that is the term we want to bound. 

We follow the ideas of Q9[] : we multiply the scheme on N d!3ab by At(N^ +l - N®) and the scheme on P d!3bb by 
At (P"^ 1 - P^) and we sum over K e T and n. It yields 



Ei + E 2 + E 3 + Fi + F 2 + F 3 = 0, 



(35) 



with 



Nt-1 



£ ' =Z Z m ^^ +1 - ^wr 1 - £2 = -Z Af Z t£ dn £> 

n=0 KeT n=0 creS 

cte£ x 

£ 3 = £ Af Z F > = Z Z m ^« +I - ^xn +1 - p ^ 

2Vr-l 

F 3 = j> £ 



n=0 o-efi 
o-e£ A ' 

2Vr- 1 



n=0 #e7~ 



n=0 cre£ 



n=0 o-efi 



As (N K +l - N" K )(N" K +l - N%) = ((N K +l - N®) 2 + (N K +l - N" K f - (N" K - N°) 2 )/2, we get: 

1 \-i n n 9 m(Q)(M-m) 2 
£, > -- J] m(X)(JV° - ^) 2 > '- 

KeT 

1 v n n i m(Q)(M - m) 2 
and F, > - - 2] ™(K)(< - P D K f > '- . 



(36) 



KeT 



We may also bound the terms E 3 and F 3 . Indeed, using successively the property of the flux T' K ^ ( IA.4ab . the L a 
estimates and Cauchy-Schwarz inequality, we get 



Nt-1 



\E 3 \ < Z Af Z ^m a x(N K + \N K ^)DAl°gN -Vy +1 D tT N D \. 



n=0 creS 
<reS K 



Nt-1 



Z A 'Z TrwmWFKNZbfaQogN-Wy" 1 ) 



n=0 <re& 
cre£ K 



1/2 



But, the right-hand-side is bounded thanks to the control of the entropy production (12 lbb and the hypothesis (0b). 
Following similar computations for F 3 , we get 



\E 3 \ < 9C(l + A 2 ) and \F 3 \ < 9C(l + A 2 ) 
with 'K depending only on T, K E , M, m, N D ', P D , /? and £ 



(37) 
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We focus now on the main terms E 2 and F%. Using the definition of the Bernoulli function (15[ , the numerical fluxes 
and ff£i, defined by |Q3), rewrite: 



Kat 



Ik 

2 



D ^(N n K +1 + Ng) - DW"£ coth 



n+l \ 



DN\ 



K.a 



s->n+ 1 _ 

^ ~ 2 
Since xcoth(x) > |x| for all x e R, we obtain 



d ^"kM P "k +1 + O - D ^ coth 



n+l \ 



DP 



K.a 



E 2 > 

n=0 



.N T -\ 

F 2 > -^Af 

n=0 



J] t„. D^" +1 (ZW" +1 ) 2 -£ T a DV£((N%f - (N" K +l ) 2 ) 



o-efi 



creE 
creE K 



Z TV ZVF" +1 (Z)^ 1 ) 2 +2 ^^((P^) 2 - (P^ +1 ) 2 ) 



cte£ 



creB 



Summing these two inequalities, we can integrate by parts due to the quasi-neutrality of the boundary conditions ([18} 
and we get 

Ni 1 . 

e 2 + p 2 > x Z Af Z Z ^^(W) 2 - (n +1 ) 2 )+ 5 r w . 



n=0 ATeT <reS K 



In the case /I = 0, using = NV' , we obtain 



E 2 + F 2 > 2^ BV ' 



(38) 



In the case /I > 0, using the scheme on *P (113c) . we get: 

1 Afr_1 1 
£2 + F 2 > — £ Af J] mCWr 1 " n +I )(W) 2 - (PT) 2 )+^Tbv. 



n=0 KeT 



Since the function x t-» x 2 is nondecreasing on R + , it also yields Q8) , Finally, we deduce the weak-BV inequality 
d22a) from (05} , (05), (07) and (08). 

4.2. Discrete L 2 (Q, T; H l )-estimates on the densities 

Now, we give the proof of the inequality (122bb of Theorem 1 1.5 1 Therefore, we start as in the proof of (122 a) with (|35V 
But, we treat in a different manner the terms E 2 and Indeed, for all K e 7~ and all cr e the Scharfetter-Gummel 
fluxes and defined by <Q3) rewrite 



^ = r^OT^)^ - K-DV&Fg-DPgyffg - T a DP% 
with P defined by B(x) = B(x) - 1 for all x e R. Therefore 

N T -1 N t -1 

E 2 + F 2 = J] Af J] T, r (ZW" +I ) 2 + Z AiJV^D^ 1 ) 2 + £ 2 + F 2 , 



(39a) 
(39b) 



(40) 



n=0 <re£ 



n=0 cre£ 



N T -1 



N T -1 



with F 2 = - Z ^ Z TZtDN'C and F 2 = - Z Af Z 



n+1 
AT.cr- 



n=0 <te8 



n=0 cre£ 
cre& K 
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But, as for the fluxes , we can rewrite the fluxes f^l either under the form (IA.2at or (I A. 2b) with B instead of B. 
Then, as x(x — y) = \{x - y) 2 + \(x 2 — y 2 ), we get either 

—^{D^f + ((N" K +1 ) 2 - (AO 2 ) 

+B(D*¥ n £)(D (r N n+1 ) 2 j, (41a) 

^{D a N n+l ) 2 - ({N n £f - (N'^f) 

+B(-D^)(D (T Af" +1 ) 2 J. (41b) 

But, B(x) > for all x < and B(-x) > for all x > 0. Then, using (l4Tal when OT^ < and (14TB when 
DT"l +1 > 0, we obtain in both cases 

-T^DN n £ > ^ (-Z^r'^ZW"^) 2 + DT"^ ((^ +1 ) 2 - (^) 2 )) . 

Similarly, we have 

-@t* Dpn £ > y (-D^" +l (D a P" +l ) 2 - DV£ ((P" K +l ) 2 - (Ptlf)) ■ 
It yields, after a discrete integration by parts 

£ 2 + F 2 > ~-T BV + 2 Z Af Z Z T - D ^ ( (A ^' )2 " (P z 1)2 ) ' 

«=() ATer o-efi A - 

and, thanks to the scheme (1 1 3cb . 

E 2 +F 2 >--T BV . (42) 

Then, we deduce the discrete L 2 (0, T; H 1 ) estimate on N and P d22bl from (J35j, ([36]l, (S3, 63, dUJi and the weak-BV 
inequality ( I22al >. 

4.5. Discrete L 2 (0, T; H x {Qj) estimate on *P 

We conclude the proof of Theorem 1 1.5 1 with the proof of the L 2 (0, 7\ //') estimate on *P ( I22cb . Once more, we use 
Theorem ll.4l in the proof. 

Let us first consider the case A — 0. In this case, multiplying the scheme on *F ( 119bt by Af( v F^. +1 - WZ) and summing 
over K eT and n e {0, . . . , N T - 1 ), we get: 

Nt-1 

£ At J] t„D¥£ (n¥& - (a^ 1 + iV#) = 0. 

n=0 <reS 

Then, thanks to the //"-estimate ( f20b and the inequality a(a — b)> a 2 /2 - b 2 /2 (Va, £> e K), we obtain: 

N T -1 N t -1 

2 A^r^CD^ 1 ) 2 < Af^r^/VF ) 2 , 

which yields d22cl ). 
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Now, let us consider the case A > 0. We follow the ideas developed by I. Gasser in I2411 at the continuous level and 
adapt them to the case of mixed boundary conditions. We set: 



N T - 1 

n=0 KeT 
N T -l 



+ Z At Z ^mm(N'^\N" K ^) + mm(Pl + \P^))(D^ n+l )\ 



n=0 o-efi 

CTE6, 



Multipliying the scheme on *P ( I13cb by Af(P^, +1 - N% + )/A~ and summing over A" e 7" and nejO,...,^-!), we get 



JVr-l 



Yj ai Tj m{K) 

n=0 ATeT 



n+1 p«+K2 ATj — 1 



(^f - p n K y 

A 2 



n=0 o-efi 
o-efi,- 



due to the quasi-neutrality of the boundary conditions (IT8V Therefore, may be splitted into J' — J7i + J2 with 

-min(^ +1 ,^)D(log^V-*P)j:A 

N T -\ , ' 

£f* 1 -(DAT^-rnin(^ +1 ,iV^)D0ogA0^) 

Applying Young inequality on Ji , we get 



N T -l 



l^ll<^ Af 



n=0 



^(Ar¥ n+ y(rrdn(A^A^) + rrdn(P£ +1 ,P^))+F +1 

< 5 Z Ar Z T -( Z V r+1 ) 2 (min(^ +1 ,^) + min^P^)) + Ke(1+ A '\ 



o-efi 
crefis: 



n=0 o-efi 
o-efit- 



Now, we estimate the term J2 which does not appear at the continuous level because VA = AV log N. For all x, y > 
we have 

, , y~ x ^ (x-y) 2 

log y - log x ■ 



min(x, y) 



2 min(x, y) 2 



It yields 



11 . 



and 



\DP'^ - mm(P" K +l , F^DQog P)^| < 



\DN n £-T^{N£\N£^D{\o E N) n £\ * 



N T -\ 



2m 



2m 



™ *± £ Ar J] r.|D^l((00 2 + (^) 2 )< ^J"^ 



n=0 o-efi;,, 

o-=a:|l 



As J' — J] + J2, the estimates on Ji and J2 imply that 

J < ^±^1(1+ ^). 
2m 

As A and P are lower bounded by m d20l >, it yields (I22cb in the case A > 0. 

16 



5. Numerical experiments 

In this section we present some numerical results in one and two space dimensions. Our purpose is to illustrate the 
stability of the fully implicit Scharfetter-Gummel scheme for all nonnegative values of the rescaled Debye length A. 

5.1. Test case 1: ID with C — 

First, we consider a one dimensional test case on Q = (0, 1), with a zero doping profile since this situation corresponds 
to the one studied in this paper. Initial data are constant No(x) = Pq(x) = 0.5, Vx e (0, 1). We consider quasi-neutral 
Dirichlet boundary conditions N D (0) = P D (0) = 0.1, ¥ D (0) = and N D (1) = P D (l) = 0.9, ^(1) = 4. 




Figure 5.1: Test case 1 . Errors in L 1 norm as functions of At, for different values of A 2 . 
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Figure 5.2: Test case 1. EiTors in L 1 norm as functions of A 2 , for different values of Ax . 

Since the exact solution of this problem is not available, we compute a reference solution on a uniform mesh made 
of 10240 = 20 x 2 9 cells, with time step At = 10 -6 , for different values of A 2 in [0,1]. This reference solution is 
then used to compute the L 1 error for the variables N, P and In order to prove the asymptotic preserving behavior 
of the scheme, we compute L l errors at time T = 0.1 for different numbers of cells 9 — 20 x 2', i e {0, 8), with 
different time steps At in [10 -5 , 10~ 2 ] and various rescaled Debye length A 2 in [0, 1]. Figure IBTTI presents the L l error 
on the electron density and on the electrostatic potential as functions of At for different values of A 2 . It clearly shows 
the uniform behavior in the limit A tends to since the convergence rate is of order 1 for all variables even for small 
values of A 2 , including zero. Similar results are obtained for the hole density. 

We plot the L l errors as functions of A 2 for different values of the space step on Figure 15.21 We still observe the 
asymptotic preserving properties of the scheme in the limit A tends to zero. Moreover, the errors are independent on 
A 2 . 
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5.2. Test case 2: ID with a discontinuous doping profile 

Here, we consider a nonzero discontinuous doping profile onfi = (0, 1), which corresponds to the physically relevant 
hypothesis, but not to the framework of our study: 



C(x) = 



-0.8 for x < 0.5, 
+0.8 for x > 0.5. 



In this case, the Poisson equation ( [Tel is replaced by 



-A 1 ^ = P-N + C. 



The initial conditions are Nq(x) — (1 + C(x))/2, Pq{x) = (1 — C(x))/2 for all x e (0, 1). And, the boundary conditions 
are still quasi-neutral and of the Dirichlet type N D (0) = 0.1, P D (0) = 0.9, *F D (0) = and N D (l) = 0.9, P D (l) = 0.1, 
*F D (1) = 4. Figure 15731 presents the error in L 1 norm between the approximate solution and the reference solution 



L 1 error on N 




L error on *P 




Figure 5.3: Test case 2. Errors in L norm as functions of At, for different values of A . 

computed as previously. We observe that the convergence rate does not depend on the value of A 2 . It seems that the 
scheme is still asymptotic preserving at the quasi-neutral limit even if the doping profile C is not zero. 



5.3. Test case 3: PN-junction in 2D 



1 r° 


P-region : 

. i 




N-region 









Figure 5.4: Geometry of the PN-junction diode 

Now, we present a test case for a geometry corresponding to a PN-junction in 2D (see Figure l5~4l >. The domain Q is 
the square (0, l) 2 . The doping profile is piecewise constant, equal to 0.8 in the N-region and -0.8 in the P-region. The 
Dirichlet boundary conditions are N D = 0.9, P D = 0.1, T° = 1.1 on {y = 0}, and N D = 0.1, P D = 0.9, ¥ D = -1.1 
on {y — 1, < x < 0.25). Elsewhere we put homogeneous Neumann boundary conditions. Initial conditions are 
N (x,y) = (1 + C(x,y))/2, and P (x,y) = (1 - C(x,y))/2 . 
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The electron density profile sat time T — 1 with a mesh made of 3584 triangles and a time step Af = 10~ 2 for A 2 — 1 
and 10~ 9 are shown in Figure 1531 We observe that the scheme remains efficient even for small values of the Debye 
length and with a large time step. 




(a) Electron density N, A 2 = 1. (b) Electron density N, A 2 = 10 9 . 



Figure 5.5: Test case 3. Electron density computed at time T = 1 with a mesh of 3584 triangles and a time step At = 10 2 for A 2 = 1 and A 2 = 10 



Appendix A. Properties of the Scharfetter-Gummel numerical fluxes 

We recall that the Scharfetter-Gummel numerical fluxes f£ + J and ffj^l. defined by (TBI can be seen respectively 
as numerical approximations of j (-VN + NV¥) ■ Vk,<t and j (-VP - PV*F) • Vk,<t on the interval [f , t n+l ). At the 
continuous level, we may rewrite -VN+NW = -NVQogN-W) and -VP-PV^V = -fV(logf + v f / ). Such equalities 
cannot be kept at the discrete level. However, we can give lower and upper bounds of J^jfl and Q"^ by terms of the 
form -JV£ +1 D(logJV - and -P^ +1 D(logP + , as shown in Proposition |Appendix AA\ 

Proposition Appendix A.l. For all K e T and all cr e &k, the flux T'^ defined by (I14ab satisfies the following 
inequalities: 

ijrn+\ 

~max(N£\N£)D(logN~ < < -min(iVf \N^)D(logN- Y)^, 

if D(logN (A. la) 



- min(iVf \N%*)D(logN- < < - max(A^ +1 , A^)D(log - 

*fD(logiV-Y)^<0. (A. lb) 

Replacing *P fey — and N by P yields similar properties for the flux G"^- 

Proof. Let K €T, first, we remark that (lA.lt is trivially satisfied if cr e £^ e „ because all the terms of the inequalities 
vanish. Let cr e &K,mt U £>1 exr since the Bernoulli function B defined by (fTBT l satisfies the property (fTTT i. the flux T^t? 
defined by ( I14ab can be either rewritten 

IF- 1 = r^D^^r 1 " B (DVg) DJV#), (A.2a) 
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or T£ = T^DVgNg - B (-D^) DN&). (A.2b) 

It implies 



' K,tr - ' o 



B(D(logN)"£)DN£ 

+ (b (DQog Ny£) - B (D^) )d/V£ 
B(-Z)aogA0^)oiV^ 

+ (b {-WogNyg) - B (-DV&) )dN"£ 
But, the definition of the Bernoulli function cfT3T > also ensures that 



and^=r, 



u *K,o-^K,(T 



Therefore, we get 



and 



log v — log x 

fi(logy-logx) = — — —x, Vx,y>0. 

y-x 



7g = ^[-DOosN-T&T 1 iB(D(logW%frB(D^))DN£], 



r£ = T„[-D(log N - 4B (-DQog Nyg}* (-DVg)) DNg] 

Now, we may use the fact that B is a non increasing function on K. Assuming that the sign of D(\og N - T)^ 1 " 1 is 
known, the sign of (b (D(log N)" K + ^) - B (D'i"'^)) and [b (-D(log N) n ^) - B (-DW^)) are also known (and oppo- 
site). Distinguishing the cases DN"^ > {N" K +l < N"^) and DN^ < (N" K +l > N"^) yields inequalities dATb . □ 



Now, we give a straightforward consequence of Proposition Appendix A. 1 as a Corollary. 

Corollary Appendix A.2. For all K eT and all cr e &k, the fluxes T^ + J and GTk\ defined by ( fT4b verify: 

D(logN- <F)£j. <-Ta mm{N n K +l ,N n K + )) (DJlogN - W)" + 'f , 
ffg D ( 1 °I p + <~r<r min(P£\P"£) (D (T (log P + V)" +1 f . 
Moreover, if min(A^ +1 , A^) > and min(P^ +1 , P n ^) > 0, we also have 



| < Tq- max(N n K +l , Ng) \D„ (log N - *F) 
|^| < T a max(/* +1 , \D a (log P + ¥) 



\H+1 | 



(A.3a) 
(A.3b) 



(A.4a) 
(A.4b) 
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